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X' 

We determine the improvement coefficients b m and 6a — 6p in quenched lattice 
QCD for a range of /3-values, which is relevant for current large scale simu- 
lations. At fixed /3, the results are rather sensitive to the precise choices of 
parameters. We therefore impose improvement conditions at constant renor- 
malized parameters, and the coefficients are then obtained as smooth functions 
of <7q. Other improvement conditions yield a different functional dependence, 
but the difference between the coefficients vanishes with a rate proportional to 
the lattice spacing. We verify this theoretical expectation in a few examples 
and are therefore confident that O(a) improvement is achieved for physical 
quantities. As a byproduct of our analysis we also obtain the finite renor- 
malization constant which relates the subtracted bare quark mass to the bare 
PC AC mass. 



1 Introduction 

Lattice QCD with Wilson quarks provides an attractive framework for study- 
ing the strong interactions beyond perturbation theory. Its main shortcoming 
consists in the explicit breaking of all axial symmetries. This is generally not 
considered a fundamental problem, but it renders the renormalization of the 
theory more complicated. In addition, the lack of chiral symmetry is at the 
origin of many 0(a) lattice artifacts in physical quantities (a is the lattice spac- 
ing), which can be rather large in practice |l],||. In order to ease continuum 
extrapolations it is therefore highly desirable to accelerate the approach to the 
continuum by cancelling at least the leading cutoff effects proportional to a. 
The theoretical framework for this goes under the name of Symanzik improve- 
ment J3|,||] , and without any loss it may be restricted to on-shell quantities and 
correlation functions at physical distances ||. 

The structure of the on-shell 0(a) improved lattice action has been de- 
rived a long time ago by Sheikholeslami and Wohlert || . A detailed discussion 
of the on-shell improved theory with Nf mass degenerate quarks can be found 
in refs. ^M. There, it was also shown how chiral symmetry may be used to 
determine some of the improvement coefficients non-perturbatively. In par- 
ticular, the improved action has been determined for both N{ = (quenched 
approximation) and JVf = 2 ]^||. Concerning the improvement of composite 
fields, non-perturbative results have been obtained for quark bilinear opera- 
tors without derivatives. So far all results are for the quenched theory. How- 
ever, while the methods used to determine the coefficients ca 0H, cy |l(| and 
C T [lljfl] a ^ so a PPly f° r -^f 7^ 0' the situation is less clear for the 6-coefficients, 



which multiply cutoff effects proportional to the subtracted bare quark mass 
m q = mo — m c . 

In the quenched approximation all the 6-coefficients can be determined 
from chiral Ward identities, by considering the more general case of mass 
non-degenerate quarks |14 , 11 , 1^] . Unfortunately these methods do not easily 



generalise to the full theory, as the theory with non-degenerate quarks requires 
many more 0(a) counterterms (for the case Nf = 3 cp. p3|). Non-perturbative 



quenched results are known for by [16], for 6 m and for the combination 6a 



op [14]. Furthermore, results for all coefficients have been published in refs. [12 



[iq|. Finally we mention that all relevant improvement coefficients are known 



to one- loop order of perturbation theory |17-2C| 

In the present paper we restrict ourselves to the determination of 6 m and 
the combination 6a — 6p in the quenched approximation. These coefficients 
are needed for the determination of renormalized quark masses using either 
the bare subtracted or the bare PC AC quark masses as starting point [21,22]. 



Motivated by this application we thus extend the previous study of ref. [14] 
to larger bare couplings, and also investigate the 0(a) ambiguities in these 
coefficients. 

The paper is organized as follows. In sect. 2 we recall the PCAC rela- 
tion and its generalization to the theory with non-degenerate quarks in the 
quenched approximation. It is then shown how the improvement coefficients 
can be isolated by defining appropriate ratios of current quark masses. Evalu- 
ating these both in perturbation theory and non-perturbatively motivates our 
strategy (sect. 3). The results, together with some checks and technical details 
are given in sect. 4, and we conclude with a summary of our findings. 

2 The PCAC relation 

In this section we assume that the reader is familiar with ref. [<J • In particular 
we refer to this reference for unexplained conventions and notation. 

2.1 Mass degenerate quarks 

Our starting point is the PCAC relation 

8^% = 2mP a , (2.1) 

where ^ is a doublet of light quarks and 

A«=^ 75 ^V, P a = ^75^V>, (2-2) 

are the isovector axial current and density respectively. On the lattice with 
Wilson quarks this relation only holds up to lattice artifacts of 0(a). These 
may be reduced to 0(a 2 ) by tuning the coefficient c sw of the Sheikholeslami- 
Wohlert term in the action, and the coefficient c\ in the improved axial cur- 
rent, 

(Ajr^A^ + CAad^, (2.3) 

such that the PCAC relation with the improved current holds exactly in a 
few selected matrix elements. In this way, non-perturbative results have been 
obtained for ca and c sw in the quenched approximation || and for c sw also in 
the case N t = 2 @. 

Given the improved action and the improved axial current, the PCAC 
relation leads to the definition of a renormalized 0(a) improved quark mass, 

z A (io A8m ) m 

Zp(l + bpam q ) 



Here, m is a bare current quark mass denned through some matrix element 
of the PCAC relation. On the other hand, the renormalized quark can be 
written in the form 

m R = Z m fh q , fh q = m q (l + b m am q ), (2.5) 

where m q = mo — m c is the subtracted bare quark mass. Combining these 
equations and expanding in powers of a we obtain 



m = ZmAl + [6 m - 6a + bp] am q ) + 0(a ), (2.6) 

where Z is a finite combination of renormalization constants, 

Z = %^, (2.7) 

which is a function of g^ = g^{\ + bgam q ) only. From this equation it is 
clear that the combination 6 m — 6a + 6p of improvement coefficients can be 
determined by studying the dependence of m upon m q at fixed go . If instead 
go is kept constant, one obtains an additional contribution, 



m = Zm q I 1 + 



2 dlnZ u 
6 m - 6a + 6 P + g 2 6 g 



am q +0(a 2 ), (2i 



where the renormalization constant is now evaluated at <7q. We note however, 
that knowledge of the combination 6 m — 6a + 6p is not sufficient to determine 
a renormalized improved quark mass, if one starts from the (measurable) sub- 
tracted bare quark mass m q = too — m c or the bare PCAC mass to. Beside 
the determination of the renormalization constants it is then necessary to de- 
termine 6 m and the combination 6a — 6p separately. 

2.2 Non-degenerate quarks and the quenched approximation 



It has been noticed in ref. 14] and later in ref. (O] that the quenched (or 
valence) approximation leads to major simplifications, which may be used to 
determine the coefficients 6a, 6p and 6 m separately. First of all the coefficient 
6 g vanishes in this approximation. But more importantly, the structure of the 
0(a) improved theory with non-degenerate quarks remains relatively simple, 
in contrast to the full theory, where non-degenerate quarks lead to a prolifer- 
ation of new improvement coefficients [O]. In the remainder of the paper we 
will restrict ourselves to the quenched approximation and shortly review the 
features that will be needed in the sequel. 



In the quenched approximation the bare quark masses are separately im- 
proved for each quark flavour, i.e. an equation of the form 

m q = m q (l + 6 m am q ), (2.9) 

holds for each flavour. The isospin symmetry is broken in the presence of 
non-degenerate masses, and it is useful to introduce the off-diagonal bilinear 
fields through 

A± = 4±m, P ± = P l ±iP 2 , (2.10) 

which are expected to satisfy the PCAC relation 

d fl A± = (m u + m d )P ± , (2.11) 

up to cutoff effects^]. In the quenched approximation the improvement of these 
off-diagonal fields is the same as in the degenerate case, except that 6 a and ftp 
now multiply the average of the subtracted bare quark masses. For instance 
the improved axial current is given by 

(A R )± = Z A [1 + 6 A ± {am^ u + am^)} (A/) J, (2.12) 

and the axial density has an analogous structure. The mass dependence is 
such that b m can be disentangled from b\ — bp by considering the analogue of 
eq. fl2.6|) for the case of non-degenerate quarks. 



2.3 Current quark masses and estimators for b m , 6a — &p and Z 

To define current quark masses we use correlation functions derived from the 



QCD Schrodinger functional [23,2 



f/(x ) = -l(A+(x)0-), / P J (*o) = -i(P + (*)0->, (2.13) 

with the source ± = O 1 ± iO 2 and 

O a = a 6 ]TC(yh5^C(z). (2.14) 



y-z 



The indices i,j refer to choices mo, u = 1^0,1 and mo,d = m o,j f° r the quark 
masses taken from a list of numerical values {mo,i,mo,2, . . . } to be specified 
later. This definition is such that the standard correlation functions /a and 



1 The flavours in the isospin doublet are generically denoted by u and d. Of course, one 
may also identify one of the flavours with the strange quark. 



/p IHIZt are recovered in the degenerate case i = j. We now define bare current 
quark masses through 

= do/A (go) + acAd$d fg (x ) 

2fAx ) 

Apart from xq these quark masses depend on the lattice size L/a, the ratio 
T/L, the angle 9 and also the precise definition of the derivatives. In 
eq. (|2.15|) we have followed ref. || by setting do = ~ (do + <9q ) , with the usual 
forward and backward derivatives given by 

dof(t) = -[/(< + a) -/(<)], (2.16) 

a 

9*of(t) = -lf(t)-f(t-a)}. (2.17) 

a 

In addition we consider current quark masses with improved derivatives [14j, 
by replacing in eq. ( |2.15 ), 



do -» d (1 - \a 2 dldo) , (2.18) 

« - % (1 - ±a 2 d* d ) . (2.19) 

When acting on smooth functions these improved lattice derivatives have er- 
rors of 0(o 4 ) only. The generalization of eq. ( |2.6| ) now reads 

rriij = Z \{m^i + m q j) + ^b^am 2 ^ + am 2 ^) 

-\(b A - b P )a(m q ,i + m qJ ) 2 ] + 0(a 2 ). (2.20) 

To isolate the coefficients we consider the combination 

2ami2 - amii - am 2 2 = /(om^ijOm^), (2-21) 

which is an analytic function of the subtracted bare quark masses. Further- 
more, it is symmetric under exchange of the arguments, f(x,y) = f(y,x), and 
vanishes for x = y, so that its series expansion can be cast in the form 

<x> 

S(x, y) = {x- y) 2 Y, cnk(x - y) 2n (x + y)\ (2.22) 

n,fc=0 

with real coefficients c n k- In particular, we note that [O] 

coo = Zl(b A -b P ), (2.23) 



up to terms of 0(a) which do not depend on the quark masses. To isolate 
6a — 6p we also consider the difference 

amu - am 22 = a(am q ,i, am q! 2), (2.24) 

and obtain an expansion of the form 



oo 
n,k=0 



ix, y) = {x-y)Y j d nk (x - y) 2n {x + y) k , (2.25) 



with doo = Z ( U P to quark mass independent terms of order a 2 ). Hence it is 
clear that the ratio 

R AP = , 2 ( 2 ""-™ii-m M ) (226) 

(mn - m 2 2){am qA - am^ 2 ) 

has a chiral limit and provides an estimate for the coefficient 6a — 6p, up to 
terms of 0(om qi i + am q , 2 ) and other quark mass independent lattice artifacts 
of 0(a). 

To obtain a similar estimate for 6 m it is useful to introduce a third quark 
mass given by the mean of the first two, 

m 0j 3 = 1(^0,1 + "10,2). (2.27) 

Taking the same steps as above we can estimate 6 m through 

= 4 <""3> (2 . 28) 

(mn - m 22 )(am q ,i - am qy2 ) 

where the error is again 0(am q i + am q ^ 2 ). We note in passing that the com- 
bination 6a — 6p — 6 m can be estimated from the ratio 

(mn - m 22 )(am q: i - am q:2 ) 

which involves only correlation functions with mass degenerate quarks (cf. 
subsect. 2.1). Furthermore one could estimate the finite renormalization con- 
stant Z through 

R z = h (6a - 6 P - 6 m )(amn + am 22 ), (2.30) 

m q ,i - m qi2 

which is correct up to 0(a 2 ) for the correct choice of 6a — 6p — 6 m . However, 
knowledge of this combination is not required as one may instead use the 
estimate ( |2.29 ). Alternative estimates of Z can be obtained from current 



quark masses which derive from correlation functions with non-degenerate 
quarks. Also in this case it is possible to cancel errors of O(am) by using the 
estimates -Rap and R m . Finally we emphasize that all the ratios have a chiral 
limit, and none of them requires knowledge of the critical mass m c . 



L/a 


o(0) 
it AP 


o(l) 

-ft AP 


o(0) 


B (l) 


o(0) 
K z 


p (i) 
K z 


standard derivatives 


8 


-0.0591 


-0.0544 


-0.4628 


-0.0819 


1.0022 


0.0881 


12 


-0.0405 


-0.0381 


-0.4753 


-0.0867 


1.0010 


0.0895 


16 


-0.0307 


-0.0297 


-0.4815 


-0.0891 


1.0005 


0.0900 


20 


-0.0248 


-0.0245 


-0.4853 


-0.0905 


1.0003 


0.0902 


24 


-0.0207 


-0.0210 


-0.4877 


-0.0915 


1.0002 


0.0903 


improved derivatives 


8 


0.0065 


0.0035 


-0.4734 


-0.0850 


0.9973 


0.0856 


12 


0.0039 


0.0024 


-0.4824 


-0.0890 


0.9988 


0.0885 


16 


0.0028 


0.0018 


-0.4869 


-0.0910 


0.9993 


0.0894 


20 


0.0022 


0.0014 


-0.4896 


-0.0921 


0.9996 


0.0898 


24 


0.0018 


0.0011 


-0.4913 


-0.0929 


0.9997 


0.0900 


oo 


0.0 


-0.0009 


-0.5 


-0.0962 


1.0 


0.0905 



Table 1: The estimates for the improvement coefficients and the renormalization constant 
in perturbation theory on lattices of size T/a x (L/a) 3 with T = 3L/2. We have set xo = T/2, 
— 0.5 and chosen the bare quark masses such that Lm^i = 0.24 and Z/m q ,2 = 0.40 for 
all lattice sizes. The expected limiting values for the coefficients are given in the last row. 
Corrections at finite lattice size are 0(a/L) for the improvement coefficients and 0(a 2 /L 2 ) 
for the renormalization constant (up to logarithms). 



2.4 Perturbation theory 

The perturbative expansion of /a and /p on a lattice of fixed size L/a has 
been explained in ref. |lq,2Cl. It is then straightforward to obtain the corre- 
sponding expansion of the current quark masses and hence of the ratios i?x 
(X = AP,m,Z), viz. 



? (0) 



2p(l) 



Rx = R'x' + 9oRx + 0(5, 



(2.31) 



After extrapolation to the chiral limit the coefficients are still affected by quark 
mass independent lattice artifacts. However, as L is then the only scale in the 
problem, these effects may be cancelled by extrapolating the lattice size L/a 
to infinity. Alternatively one may fix the subtracted bare quark masses in 
units of L, and then take the limit L/a — > oo, i.e. the chiral and the infinite 
volume limit are reached at the same time. An example is given in table 1. 
As expected the estimates for the coefficients converge to the known results 
of ref. p0|| , which were obtained by requiring scaling of renormalized finite 
volume correlation functions. For gauge group SU(3) and neglecting terms of 



0(<7q) one has 

b m = -\ - 0.09623(3) x g%, (2.32) 

b A - 6 P = -0.00093(8) x g\. (2.33) 

The perturbative result for the renormalization constant Z is implicit in the 
literature |jj| and has been computed by one of the authors |2(| . For N = 3 
it is given by 

Z = 1 + 0.090514(2) x gl + O(^). (2.34) 

3 Improvement conditions 

3.1 Improvement conditions at fixed lattice size 

Beyond perturbation theory the estimates for the improvement coefficients are 
affected by an intrinsic ambiguity of order a/r^ (tq denotes the hadronic scale 
of ref. [§7]]). In contrast to the perturbative example of the last section it 
is thus impossible to eliminate this ambiguity by chiral and infinite volume 
extrapolations. Therefore the estimates Rx at some fixed lattice size L/a, at 
fixed quark masses amo^ and for a definite choice of T/L and 6 may be taken as 
a definition of the improvement coefficients, provided the 0(om q ) and 0(a/L) 
lattice artifacts are not too large. At our strongest coupling j3 = §jg\ = 6.0 
we expect typical O(a/ro) ambiguities to be around 0.1 or less. Taking this 
as a guideline we use perturbation theory to choose the lattice size and the 
quark mass parameters. 

A typical example is summarized in table 2. The lattice size is 12 x 8 3 , 
the angle 9 = 0.5, and the bare masses are om qi i = 0.03 and am qt 2 = 0.05. To 
obtain the one-loop coefficients we have used c A = —0.00757 and the critical 



quark mass has been set to its value in infinite volume [18]. However, as 
expected, the dependence upon the latter is weak. We note that the improved 
derivative yields consistently better results also at the one-loop level, where 
this is a priori not expected. With the chosen parameters, the perturbative 
results show the expected behaviour for the cutoff effects. In particular, we 
have checked that the deviation from the exact results becomes smaller with 
decreasing quark masses. 

Evaluating the same estimators non-perturbatively at (3 = 6.0, with ap- 
proximately the same choice of parameters (cp. sect. 4), and with the non- 
perturbative values for c sw and ca ||, we obtain the results shown in figs. || 
and y. Somewhat surprisingly the non-perturbative results for b^ — bp are quite 
sensitive to the choice of parameters, and to the choice of the lattice derivative. 



Xq 


o(0) 


o(l) 

-ft Ap 


o(0) 
-n-m 


p(i) 

-Km 


o(0) 
K z 


o(l) 

K z 


standard derivatives 


3 


-0.052 


-0.035 


-0.456 


-0.086 


1.002 


0.089 


4 


-0.054 


-0.039 


-0.458 


-0.082 


1.002 


0.090 


5 


-0.056 


-0.046 


-0.460 


-0.081 


1.002 


0.089 


6 


-0.059 


-0.054 


-0.463 


-0.082 


1.002 


0.088 


7 


-0.062 


-0.062 


-0.466 


-0.084 


1.002 


0.087 


8 


-0.065 


-0.068 


-0.469 


-0.085 


1.002 


0.085 


9 


-0.069 


-0.069 


-0.472 


-0.086 


1.002 


0.084 


improved derivatives 


3 


0.009 


0.020 


-0.471 


-0.098 


0.997 


0.094 


4 


0.008 


0.005 


-0.472 


-0.087 


0.997 


0.088 


5 


0.007 


0.004 


-0.473 


-0.086 


0.997 


0.086 


6 


0.007 


0.004 


-0.473 


-0.085 


0.997 


0.086 


7 


0.006 


0.003 


-0.474 


-0.085 


0.997 


0.085 


8 


0.005 


0.002 


-0.475 


-0.085 


0.997 


0.085 


9 


0.003 


0.011 


-0.477 


-0.086 


0.997 


0.089 



Table 2: The estimates for the improvement coefficients and the renormalization constant 
Z in perturbation theory on a f 2 x 8 3 lattice. The parameters are as explained in the text. 
The one-loop coefficients have been evaluated for gauge group SU(3). 



The situation may be slightly improved by correcting for the tree-level cutoff 
effects, i.e. by setting 



b A -b P = R A p - R 



(0) 
AP' 



Rn 



R® 



i 

2' 



(3.1) 



such that the estimators assume the correct values in the continuum limit. 
However, from the perturbative results of table 2 we infer that this only cor- 
rects for a fraction of the discrepancy in the case of 6a — bp, while the xq 
dependence of the estimator for 6 m becomes somewhat more pronounced. 



3.2 Improvement conditions at constant physics 

At least the case of 6 a — bp seems different from the determination of c sw and 
ca in ref. ||, where the results were quite stable against a variation of the 
parameters. Already in perturbation theory the O(a) effects in 6a — 6p are 
much larger than in c sw and ca- One might argue that our quark masses are 
not small enough, or that the lattice size is too small. However, it is our general 
experience that the spread of the estimates is considerably larger than naively 
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Figure 1: Estimates of &a — bp versus Xo- Squares are obtained with the standard, circles 
with the improved derivative. 



expected, also on larger lattices or if a chiral extrapolation is attempted. It 
should be emphasized at this point that this does not imply a failure of the 
improvement programme. Improvement is an asymptotic concept, and thus 
only determines the rate of the continuum approach. It is a priori not clear 
how large the intrinsic 0(a) ambiguities typically are for a given improvement 
coefficient. 

There are essentially two alternative ways of dealing with this situation. 
First one may determine a "typical" spread of the values for Rx (X ^ Z) 
and quote a corresponding systematic error for the estimate of the improve- 
ment coefficient. In this paper we follow a second approach by imposing the 
improvement condition at constant physics. More precisely, we fix all renor- 
malized quantities in units of L and keep L/tq constant as we approach the 



10 



-0.6 



-i 1 1 1 r 



-i 1 1 1 r 



-0.65 



i * 



.a -0.7 



-0.75 



-0.8 



_i I i i i L 



_i i i i i_ 



6 



10 



Figure 2: Estimates of 6 m versus Xo- Squares are obtained with the standard, circles with 
the improved derivative. 



continuum limit. As a consequence the estimates become smooth functions of 
a/ro and thus also of go- Furthermore it is obvious that any other estimate 
Rx may have a different dependence upon go, but the difference Rx — Rx is 
again a smooth function which must vanish in the continuum limit with a rate 
proportional to a/ro. Note that the very same strategy has previously been 
applied to the determination of the finite renormalization constants for the 



isospin currents [16]. However, in the case of renormalization constants the 
intrinsic ambiguity is of 0(a 2 ). 

Our procedure highlights the importance of the continuum limit. Fur- 
thermore we emphasize that in cases where the intrinsic O(a) ambiguity is 
large, the important result is not so much a numerical value for the improve- 
ment coefficient at fixed (3 but rather its dependence on go which follows from 



11 



keeping the physics fixed. 
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Figure 3: This plot illustrates how well the physics is kept constant. Full symbols refer to 
the bare current quark masses, empty symbols to the renormalized quark masses (the latter 
are not available at the largest (5 value). 



4 Technical details and results 

We describe the technical aspects of the numerical simulations in detail and 
present our results. 

4.1 Choice of simulation parameters 

The improved action and axial current have been determined for f3 = 6/<7q > 
6.0 B. This sets the upper limit for the bare couplings to be considered. At 



12 



(3 = 6.0 the probability of encountering a quark zero mode (so-called "ex- 
ceptional configurations") increases rapidly if the lattice volume is too large 
or the quark mass too small ||. To avoid this problem we choose a 12 x 8 3 
lattice, 9 = 0.5 and quark masses such that am\\ = 0.03 and am,22 = 0.05. 
The parameters are thus approximately the same as in the perturbative exam- 
ple discussed in the preceding section. Note that by using the current quark 
masses instead of the subtracted bare quark masses we avoid the determina- 
tion of the critical quark mass. In fact, as none of our observables requires the 
knowledge of m c we will not determine this parameter at all. 

We want to keep the physics constant as we vary (3 towards larger values. 
To this end we fix all dimensionfull quantities in units of the physical scale 
tq [27 1, which has been determined very precisely in ref. |28(| . From the fit 
function, 

ln(a/r ) = -1.6805 - 1.7139(/3 - 6) + 0.8155(/3 - 6) 2 

- 0.6667(/3 - 6) 3 , (4.1) 

which is valid in the interval 5.7 < j3 < 6.57, we infer that our lattice size L/a = 
8 at (5 = 6.0 corresponds to L/vq = 1.49. The requirement that this ratio be 
constant as the lattice size L/a is increased determines the corresponding (5 
values. Furthermore, we choose T = 3L/2 and tune the bare quark masses 
mo,i and too,2 such as to maintain 

Lm n = 0.24, Lm 22 = 0.40 (4.2) 



for all lattice sizes with the third bare quark mass then given by ( 2.27| ). The 



resulting simulation parameters are summarized in table |3[ Note that we 
quote the values for the hopping parameter k = l/(2amo + 8) rather than 
amo, since we use this parameter in our programs. A complication arises at 
our largest lattice size, L/a = 24, as the corresponding (3 value lies clearly 
outside the range of validity of eq. (f4,l|) . In order to include the largest lattice 
we have extrapolated the data points in J^] using a linear fit of ln(a/ro) as 
a function of (3. Of course this procedure is not unique. From the spread of 
results obtained with different ansatze for the fit we estimate the error in (5 to 
be at most 0.003. However, from this procedure it is not possible to estimate 
the real error in a reliable way and the results obtained at (3 = 6.756 have to be 
taken with some care. In fig. |3| one can see how precisely the conditions ( [4.2[) 
are satisfied. The relative accuracy is better than 2 per cent. The reader may 
be worried that we are keeping fixed the bare current quark mass in units 
of L, rather than corresponding renormalized masses. It turns out that the 
difference is small in the range of (5 values considered. To illustrate this, we 
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L/a 


(3 


«l 


K2 


8 


6.0 


0.133931 


0.133241 


10 


6.13826 


0.134828 


0.134263 


12 


6.26229 


0.135107 


0.134654 


16 


6.46895 


0.135126 


0.134794 


24 


6.756 


0.134856 


0.134634 



Table 3: Values for f3, ki and K2 for simulations with constant physical volume and with 
Lmu = 0.24 and Lm 2 2 = 0.4. 



have, for L/a < 16, also plotted Lm R = LmZj^/Zp for both masses. They 



are obtained using the non-perturbative result for Z\ [16] and Zp in the SF 



scheme [29,21]. The largest lattice had to be left out, as a non-perturbative 
result for Zp is not available for (3 = 6.756. As the renormalization constant 
barely varies over the range of couplings considered, also the renormalized 
quark masses are constant to a good precision. 



4.2 The numerical simulation 

Our numerical simulations were performed on APE/Quadrics parallel comput- 
ers with SIMD architecture and single precision arithmetic following the IEEE 
standard. We have used machines with up to 512 processors with an approx- 
imate peak performance of 50 MFlops per node. While we have distributed 
the larger lattices over the whole machine the smaller lattices have been simu- 
lated with up to 32 independent replica. To generate the gauge fields we have 
applied a hybrid overrelaxation algorithm with two overrelaxation steps per 
heatbath sweep. Every 15th complete update step the fermionic correlation 
functions have been measured for all combinations of quark masses. With 
these choices the integrated autocorrelation time between successive measure- 
ments of the fermionic correlations was found to be negligible. All in all we 
have accumulated between 870 and 3712 independent measurements for the 
various lattice sizes. From expectation values of the fermionic correlation func- 
tions our secondary observables (|2.26| ), (|2.28| ) and (|2.30| ) can be constructed 
as described in sect. 2. To increase statistics we have chosen to average our 
secondary observables over the time slices L/(2a), . . . , (T — L/2)/a. Note that 
the number of time slices used for the average is scaled with the lattice size, 
so that also in this respect the requirement of constant physics is satisfied. To 
compute the errors of b\ — bp, b m and Z we have calculated the autocorrelation 
functions of these observables along the lines of appendix A in [3C]. 
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p 


9l 


b A -b P 


b m 


Z 


6.0 


1.0 


0.1703(49) 


-0.7127(48)(7) 


1.0605(4) 


6.13826 


0.9775 


0.0563(23) 


-0.6915(35)(11) 


1.0907(4) 


6.26229 


0.9581 


0.0251(33) 


-0.6852(36)(30) 


1.0992(4) 


6.46895 


0.9275 


0.0090(52) 


-0.6846(44) (68) 


1.1048(3) 


6.756 


0.8881 


0.0007(36) 




1.1047(2) 



Table 4: Non perturbative results for b\ — bp, b m and Z. We have averaged over the 
time slices L/(2a), ... , (T — L/2)/a, and quote the corresponding statistical errors. In the 
case of 6 m , the second error is an estimate of the uncertainty due to rounding errors in the 
determination of the third n value. Further details are given in the subsect. 4.4. 



4.3 Result for 6a — bp 

Our results are tabulated in the second column of table HI and shown in 
figure |3|. At small values of the bare coupling the non perturbative result is 
consistent with the one-loop result ( |2.33D . At larger couplings 6a — bp increases 
rapidly. Our numerical results are well described by the function 



(bA-b P )(g 2 ) 



2 1 + 23.3060 gk - 27.3712 gft 
-0.00093 9 2 x '" "' 



1 - 0.9833 g$ 



(4.3) 



which incorporates the perturbative result for small values of g^. In the range 
0.8881 < 5q — 1 our data are represented with an absolute deviation smaller 
than 0.003. 

4.4 Result for 6 m 

The determination of the improvement coefficient 6 m requires a subtle can- 
cellation which is achieved by setting the third bare mass parameter to the 



average of the other two ( 2.27 ). In terms of the hopping parameters this 
translates to 



«3 



2«i«2 
«1 + K 2 



(4.4) 



As we perform simulations with single precision arithmetic following the IEEE 
standard we expect relative roundoff errors of the size 6x 10 -8 . This introduces 
a systematic error which however can be estimated by repeating the derivation 
of formula fl2.28|) and taking into account a small deviation of |(mo,i +^-0,2) — 
7710,3 from zero. This difference can be parametrized by the corresponding 
hopping parameters which are the input parameters of the simulation program. 
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Figure 4: Results for 6a — bp as a function of g$ from numerical simulations (filled squares) 
and one-loop perturbation theory (dashed line). The full line represents the fit (4.3). 
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Thus we define 






K 3 . 



(4.5) 



Then an additional term proportional to e appears in ( 2.28; ). The result is 

\2 



R 



corrected 



Rm + 



2e 



K1 + K 2 ) 



(4.6) 



(«i - k 2 ) 2 ' 

up to terms of 0(e 2 ). This effect is tabulated in table [|. To obtain these values 
we have used the single precision internal values for the hopping parameters. 
The correction constant turns out to be negative for all lattice sizes except 



L/a 


8 


10 


12 


16 


24 


Sbm 


-0.0013 


-0.0021 


-0.0059 


-0.0135 


0.040 



Table 5: Correction constant for 6 m as computed from equation (4.6) 



for L/a = 24. While on the smaller lattices it is rather small it grows with 
the lattice size since the masses and thus k\ and k 2 get closer to each other. 
We have decided to correct the statistically obtained values of b m by adding 
<5(, m for lattice sizes up to L/a = 16. At L/a = 24 this effect is about five 
times larger than the target value for the statistical error of b m so that we 
have chosen to omit this lattice in the discussion of b m . As an estimate for 
the remaining systematic uncertainty we quote 2<5ft m . Our final results for b m 
are given in table H and shown in figure 0. For large values of the coupling go 
the value for b m decreases rapidly. We find that the data are well described 
by the rational fit function 

0.6905 gl + 0.0584 g% 



M<7o) = (-0-5 - 0.09623 g 2 Q ) 



1 



1 



0.6905 gl 

which asymptotically reproduces the perturbative result ( p. 3^ 
scribes the data with an absolute deviation smaller than 0.013 



(4.7) 



The fit de- 



4.5 Alternative improvement conditions 

In order to verify the expectation that the difference between two determina- 
tions of the improvement conditions vanishes with a rate proportional to a/r§ 
we now consider a few examples. We start with b\ — bp, where the choice of 
the lattice derivative (improved vs. standard, cf. sect. 2) made a big difference 
(cf. fig. |l]). We thus define the difference 



A(6a — bp) = i?Ap(hnpr. deriv.) — i?Ap(std. deriv.) 



(4.8) 
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Figure 5: Results for b m as a function of go from numerical simulations (filled sq uare s) and 
from one- loop perturbation theory (dashed line). The full line represents the fit (|4.7|). 
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and otherwise do not change the parameters. The result is shown in fig. || 
The behaviour is approximately linear in o/ro- 

The same difference for b m is rather small, cp. fig. ^. Also alternative 
definitions of A6 m such as 



A6 m = R m (x = T/2) - R m (x = T/3). 
yield small values for the 0(a) ambiguities of b m . 



(4.9) 



0.3 - 



J" 0.2 



0.1 



-i 1 1 r 



~i 1 1 r~ 



_i i i i_ 



_L 



_i i i i_ 



0.05 



0.1 
a/r 



0.15 



0.2 



Figure 6: The difference A(6a — bp) versus a/ro- 



4.6 Determination of Z 

For the determination of the renormalization constant in the O(a) improved 
theory we would like to keep the physics fixed up to errors of 0(a 2 ). Fortu- 
nately this is already the case, for the smallness of 6a — bp implies that the 
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renormalized 0(a) improved quark masses are constant in units of L with good 
numerical precision (cp. fig. pi). To evaluate the estimator Rz one needs the 
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1.1 - 



N 1.08 - 
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1.04 
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0.9 



0.92 



0.94 



0.96 0.98 



Figure 7: Results for Z as a function of g\ from numerical simulations (filled squa res) a nd 
1-loop perturbation theory (dashed line). The solid line represents the fit function (4.1C). 



combination of improvement coefficients bp^ — b-p — b m [cp. eq. ( 2.30J) ]. We here 
simply insert the value obtained with the fit functions (4.3) and (|4.7|) . At the 
largest lattice size of our simulations go is outside the range of validity of ( [4.7D . 
We use a linear fit to the last three data points for b m to extrapolate our data 
to this point. The resulting systematic error can be estimated by calculating 
Z with different values for b m . We find that it is smaller than 10~ 4 and will 
neglect it in the following. Our results for Z are summarized in the fourth 
column of table and shown in figure p]. In the range of couplings considered 
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Figure 8: This plot shows the difference between the estimates of Z obtained from two 
alternative renormalization conditions. 



we find that our data is represented by the fit function 

1 - 0.9678 g 2 + 0.04284 g\ - 0.04373 g% 



Z{gl) = {l + 0.090514 g 2 ) 



1 - 0.9678 g$ 



(4.10) 



with a relative precision better than 0.04%. Again the perturbative result Q2.34{ ) 
has been incorporated into the fit. 

It is important to check that another determination of Z differs by 0(a 2 ). 
In fig. H we plot the difference between Rz determined with the improved and 
with the standard lattice derivatives versus a 2 jr\ and see the expected roughly 
linear dependence. 
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5 Summary and Conclusion 

We have carried out a detailed non-perturbative determination of the improve- 
ment coefficients 6 m and 6a — 6p for the range of (5 values used in present large 
scale simulations. In the case of 6a— 6p we find that the typical 0(a) ambiguity 
at j3 = 6.0 is not small. Rather than quoting a corresponding systematic error 
we advocated the use of improvement conditions at constant physics. This 
generalizes the strategy previously applied to renormalization constants [16]. 



As a result, we have obtained the improvement coefficients as smooth func- 
tions of the bare coupling [eqs. (|4.3|) , (|4.7|) 1, and we checked that other choices 
of improvement condition lead to differences which vanish with a rate propor- 
tional to a I Vo- As a byproduct of our approach we also obtained the finite 
renormalization constant Z = Z m Zp/ZA, eq. ( [4.10 ). The same remarks as 



above apply to this case, except that the difference to other determinations 
vanishes with a rate proportional to a 2 /r 2 . In contrast to the improvement 
coefficients, our result for Z differs from 1-loop perturbation theory in the bare 
coupling by no more than 3%. 

Imposing physical improvement conditions and the error analysis are fa- 
cilitated by the use of the Schrodinger functional technology and the use of 
simple direct estimators. It would be more difficult to ensure the constant 
physics condition if smearing techniques were necessary to enhance the signal, 
or if the improvement coefficients were obtained through a fitting procedure. 
Nevertheless, where they overlap, our new results are compatible with those 



of ref. [14] within the 0(a) and 0(a 2 ) ambiguities discussed in this paper. A 
drawback of the present method is that in the case of b m we reached the limit 
of single precision arithmetic, because the estimator of this coefficient relies on 
the cancellation between bare quark masses which only happens up to round- 
ing errors. As explained in detail, this problem becomes more pronounced as 
the lattice size increases and prevented us from using our largest lattice for 
an estimate of 6 m . We note that the strategy proposed here may be applied 
in other cases and may certainly be combined with the ideas of refs. [0,0] in 
order to determine the coefficients 6a and 6p separately. 

We end with some comments concerning the coefficient ca and the mag- 
nitude of the corrections associated with 6a — 6p and 6 m . Since the original 
non-perturbative determination in ref. ^ the relatively large value of c\ at 
(5 = 6.0 (as compared to one-loop perturbation theory) has been subject to 
doubts by other authors |lljl2|. In these papers, an alternative determination 
of ca was presented which yields a numerical value which is roughly half as 
big at (3 = 6.0.FI One should note, however, that all the scaling tests that have 

2 Some time ago, we had also considered alternative improvement conditions to determine 
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been carried out for physical quantities |32 : 3JJ show the expected 0(a) im- 
proved continuum approach, albeit with sizable 0(a 2 ) artifacts in some cases. 
Our present work contains additional examples of scaling. In fact, the quan- 
tities shown in figures g and || have the advantage that they constitute pure 
lattice artifacts (rather than having an a priori unknown continuum limit). 
More precisely, A6 m and A (6a — 6p) have to vanish in the limit a/ro = only 
when the theory is order a improved, i.e. when c sw and ca have the proper 
values. Similarly AZ = 0(a 2 ) is true only in this case. The convincing agree- 
ment with the expected behavior (cf. figs, y and ||) therefore constitutes new 
evidence that 0(a) improvement has been correctly implemented. 

We have emphasized above, that in taking the continuum limit of phys- 
ical quantities the correct dependence of the improvement coefficients on go 
is crucial. In addition, their overall magnitude determines whether they are 
relevant in practical applications. As we are considering coefficients multiply- 
ing quark masses, this depends in particular on whether one is interested in 
the physics of light quarks or heavy quarks. In the range of bare couplings 
considered, the bare strange quark mass in lattice units is smaller than 0.04 
|22] . Therefore, approximating the 6-coefficients determined here by 1-loop 
perturbation theory (as has been done in 122]) does not introduce errors be- 
yond the percent level for light quarks. On the other hand, for quark masses 
around the charm quark mass, non-perturbative values should be used if one is 
interested in precisions of a few percent. In addition, the difference A(6a — 6p) 
translates into an effect of about 10% for the situation of a D-meson |34| at 
(3 = 6.0 and about half as much at [3 = 6.2. Unless one takes the continuum 
limit, this is the magnitude of 0(a 2 m q ) effects to be expected for D-mesons in 
the improved theory and one may suspect 0(a 2 m 2 ) terms to be even larger. 

This work is part of the ALPHA collaboration research programme. We 
would like to thank G. de Divitiis and M. Kurth for collaboration in the early 
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edged. S. Sint acknowledges support by the European Commission under grant 
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Forschungsgemeinschaft (Graduiertenkolleg GK271). 



ca [fell]. In all cases when the lattice artifact had a large sensitivity to ca, the results were 
consistent with the original determination of ref . fel . 
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